Dynamical stabilization of matter-wave solitons revisited 
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We consider dynamical stabilization of Bose-Einstein condensates (BEC) by time- 
dependent modulation of the scattering length. The problem has been studied before by 
several methods: Gaussian variational approximation, the method of moments, method of 
modulated Townes soliton, and the direct averaging of the Gross-Pitaevskii (GP) equation. 
We summarize these methods and find that the numerically obtained stabilized solution 
has different configuration than that assumed by the theoretical methods (in particular a 
phase of the wavefunction is not quadratic with r). We show that there is presently no 
clear evidence for stabilization in a strict sense, because in the numerical experiments only 
metastable (slowly decaying) solutions have been obtained. In other words, neither numer- 
ical nor mathematical evidence for a new kind of soliton solutions have been revealed so 
far. The existence of the metastable solutions is nevertheless an interesting and complicated 
phenomenon on its own. We try some non-Gaussian variational trial functions to obtain 
better predictions for the critical nonlinearity g cr for metastabilization but other dynamical 
properties of the solutions remain difficult to predict. 



I. INTRODUCTION 

The nonlinear Schrodinger equation (NLSE) 
appears in many models of mathematical physics 
and has numerous applications. The one- 
dimensional NLSE is famous due to its in- 
tegr ability and soliton solutions. The two- 
dimensional and three-dimensional versions do 
not have such properties and are much less ex- 
plored. 

In the last decade dynamics of BECs has at- 
tracted enormous amount of interest which in 
turn is causing a renewed growth of interest in 



the NLSE, since it is known that NLSE (of- 
ten called the Gross-Pitaevskii (GP) equation in 
that context) describes the dynamics of BEC at 
zero temperature very well jjj. 

While early analytical studies of BECs were 
concentrated on (quasi-)one-dimensional sys- 
tems, (quasi-)2D and 3D systems are more im- 
portant for real experiments. In 2D and 3D sys- 
tems analytical treatment of NLSE is very diffi- 
cult and one has to use approximate methods. 

One of the very interesting and complicated 
phenomena being studied recently is stabiliza- 
tion of BEC by the oscillating scattering length 
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in two and three dimensions. eluded from the dynamics, say, due to a tight 

In ID geometry, bright solitons can be stable confinement). The system is described by the 

without trapping potential if nonlinearity is at- GP equation: 

tractive and sufficiently strong. In NLSE with q^u \ u 2 (t) 

i^ = --VV + -^rV + 0(t)MV, (I) 

attractive interaction (corresponding to BEC at z z 

with negative scattering length) in 2D free space, where r 2 = x 2 + y 2 and g(t) = 

kinetic energy can balance interaction energy at (&-Kmuj z /h) 1 l 2 Na{t) describes the strength 

certain critical value of nonlinearity g cr , but the of the two-body interaction. The interaction 

resulting solution (Townes soliton) is unstable. g(t) is rapidly oscillating: g(t) = g Q + gi sin(Qi), 

That is, if nonlinearity is either increased or de- while the confinement trap described by u r (t) 

creased (and kept fixed afterwards), the solu- is slowly turned off. Refs. J, B ^ Q suggest 

tion either expands or collapses correspondingly, it is possible to obtain a dynamically stabilized 

It was shown by several authors that stabilized bright soliton in free space in such a way. 

solutions are possible with the oscillating scat- Interactions between such objects were very 

tering length. The oscillations of the scattering recently studied in Ref. |jg. This is a very 

length lead to creation of pulsating condensate, interesting phenomenon not only in the context 

i.e. some kind of breather solution. One can of BECs but also from a broader scope of 

draw an analogy with Kapitza pendulum (a pen- nonlinear physics. 

dulum with a rapidly oscillating pivot), where Such kind of stabilization in 3D has also been 

unstable equilibria of unperturbed system is sta- reported The latter finding is, however, in 

bilized by means of fast modulation. This idea some disagreement with other investigations on 

was already applied to stabilization of beams in this topic (for example, Ref. 0). In Ref. [ill it 

nonlinear media [2j . Among many other appli- was shown that the scattering length modulation 

cations in related fields, the atom wire trap sug- may indeed provide for the stabilization in 3D, 

gested in Ref. |3| should be mentioned. In Refs. but only in combination with a quasi-lD peri- 

[J, |a] the novel application of this stabilization odic potential. So 3D geometry might need addi- 

mechanism to BEC physics was presented which tional careful examination. In the present paper 

in turn encouraged several other works on that we concentrate on quasi-2D case only, where also 

. nnnnn 

subject p, 13, \Q, llil- LLi| • not everything is clear yet. Unlike conventional 

We consider here the problem of stabilization ID solitons, higher-dimensional solitonic objects 

of BEC in 2D free space by means of rapid os- may decay. Therefore, it is interesting to inves- 

cillations of the scattering length in a greater tigate the following question: is there indeed a 

detail (the third dimension is assumed to be ex- novel genuine breather solution behind the phe- 



3 



nomenon of stabilization? As we show in this lized wavefunction does not have such parabolic 
paper, it turns out that the phenomenon does phase factor and does not fit into self-similar 
not fit into simple models being suggested ear- patterns implied by the above-mentioned meth- 
lier. For theoretical description of the process, ods. This qualitative difference between the ex- 
several methods were used by different groups of act numerical solution and all theoretical mod- 
authors: variationa approximation based on the els considered so far was not mentioned earlier. 
Gaussian anzatz [4, Q], direct averaging of the Besides, we noticed presence of steady outgoing 
GP equation a method based on modulated flux of atoms in numerical stabilized solutions. 
Townes soliton p], and the method of moments So, even numerically there is no 2D soliton so 
P). Surprisingly, we find all the methods are not far, but some slowly decaying object instead, 
very satisfactory even for qualitative predictions. Section 2 reviews the abovementioned theoret- 
In brief, the direct averaging of the GP equation ical methods. In Section 3 we give some re- 
has the disadvantage of ommiting terms which suits obtained using the variational approxima- 
are of the same order as those responsible for tion with non-gaussian trial functions, including 
creation of the effective potential, while three " supergaussian anzatz" . It is shown that a bet- 
other methods, although very different, all rely ter accuracy can be obtained for predicting crit- 
on the unwarranted assumption of parabolic de- ical nonlinearity g cr , but we were not able to de- 
pendence of the phase of the stabilized wavefunc- termine accurately such dynamical properties as 
tion on r: arg ip = a(t) + [5{t)r 2 . We find that the frequency of slow oscillations. Additionally, 
the behavior of the exact numerical wavefunc- we checked the supergaussian anzatz for another 
tion is, however, completely different (see Fig. problem: determination of critical number of at- 
2j). The above-mentioned parabolic approxima- tractive BEC in a parabolic well, and found it to 
tion (PA) of the phase factor is very popular be much more accurate than the usual gaussian 
because it is appealingly simple and indeed of- anzatz. This example also demonstrates that 
ten appears in solutions of the time-dependent the stabilization mechanism is essentially more 

n 

GP equation [13j. Usually it comes from self- complicated than that assumed by the present 

similar time evolution of the condensate den- (PA-based) methods, because predictions of the 

sity, for example in 3D the following dynamics supergaussian anzatz for dynamical properties 

of the condensate density is possible p(x,y,z) = of the stabilized solution are much less accurate 

[\i(t))\2(t)\3(t)]~ 1 p(x/ \i(t), y/\2(t), z/\s(t)) , than in static problems. 

where coefficients A; are coupled by nonlinear In Section 4 numerical results are presented 

differential equations. It is the important finding and com pared with predictions of the theoretical 

of the present paper that in our problem a stabi- me thods discussed in Sections 2 and 3. Configu- 
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ration of stabilized solution is discussed and dy- 
namics of some integral quantities of the solution 
is investigated. 

In Section 5 concluding remarks are given. 
We mention the relation between the BEC stabi- 
lization problem and stabilization of optical soli- 
tons in a layered medium with sign-alternating 
Kerr nonlinearity. 

II. SEVERAL APPROXIMATE 
METHODS TO STUDY THE PROBLEM: 

PA-BASED METHODS (GAUSSIAN 
VARIATIONAL APPROXIMATION, THE 
MODULATED TOWNES SOLITON, THE 
METHOD OF MOMENTS), AND THE 
DIRECT AVERAGING OF THE GP 
EQUATION. 

A. PA-based methods 

1. Gaussian variational approximation 

The variational approach based on the Gaus- 
sian approximation (GA) is one of the most of- 
ten used in studying dynamics of the GP equa- 
tion. In actual calculations this approximation 
however often gives a large error as compared 
to exact numerical results [7, 121. For example, 

n 

in Ref. [12|] the Gaussian approximation in dy- 
namics of attractive BEC was compared to ex- 
act numerical solution of the GP equation. It 
was found that in estimating the critical num- 
ber 7V C of the condensate (the maximal number 
of condensed particles in a trap before collapse 
occurs) the Gaussian approximation gives a 17% 



error, and similar values of discrepancy for other 
dynamical quantities (as a useful test, in the 
Appendix we provide corresponding results ob- 
tained with a supergaussian variational ansatz). 
However, it seems that in this example GA en- 
ables to reproduce important features of the sys- 
tem at least qualitatively. The GA was also used 
in many other treatments of the GP equation 
using a variational technique. In particular, it 
was applied to the problem of BEC stabiliza- 
tion by the oscillating scattering length. The La- 
grangian density corresponding to the GP equa- 
tion © is 



L\ 



dt V dt ' 



K 


dtp 




dr 



(2) 



The normalization condition for the wave- 



function is 2tt Jq 



*rdr = 1. 



In Ref. [J], a variational method with the 
following Gaussian anzatz was used, 



ip(r, t) 



1 



y/7TR(t) 



cxp 



. R(t) 2 
+ % „ „ , - r 



2R?{t) 2R(t) 



(3) 



where R(t) is the variational parameter that 
characterizes the size of the condensate, and the 
phase factor of the wavefunction describes the 

After substitution of expression © into the 
Lagrangian density © one obtains the effective 
Lagrangian L = 2ir rL[ip]dr and the corre- 
sponding Euler-Lagrange equations of motion. 
One can obtain then the equation of motion for 



5 



R(t) as dictions © and © based on the Gaussian ap- 

l g _|_ g lS \ n Qt proximation can catch (gi/Q) 1 / 2 dependence of 

R 3 (t) 2tt R s (t) the monopole moment < r > and (tl/gi) de- 

So the gist of the model is to represent the pendence of the breathing-mode frequency ujbr 

2D BEC as a classical nonlinear pendulum with but cannot determine the corresponding coeffi- 

modulated parameters. It is important that cients of proportionality, of which the one in © 

other one-parameter PA-based anzatzes also becomes infinity while the one in © becomes 

give the same nonlinear pendulum (R = (a + zero for g = -2ir, the value actually used in 

bsin£lt)/R 3 , where a,b depend on the parame- the numerical calculations. On the other hand, 

ters gi,g ,Q), but with different functional de- from numerical calculations they were able to 

pendence of a,b on the parameters. determine the coefficients as 1.06 and 0.32 cor- 

The authors of Ref. jj] use then the Kapitza respondingly (see Fit* 2 of Ref. |4j]). It was also 

averaging method to study behavior of the determined in Ref. 1.4] that in order to stabilize 

system with the rapidly oscillating scattering the bright soliton, \g \ must exceed the critical 

length. They assume the dynamics of R can be value of collapse \g cr \. Their numerical estimate 

separated into a slow part Ro and a small rapidly for \g cr \ is « 5.8 while theoretical estimate based 

oscillating component p: R = Ro(t) + p(flt). on Gaussian approximation is 2ir 6.28. The 

From the equations of motion for R and p one 2tt estimate in fact corresponds to fitting the so- 

extracts the effective potential for the slow vari- called Townes soliton by a Gaussian trial func- 

able U(Rq) ~ 4f + 4l and determines its mini- tion as will be discussed below, 

mum Inspired by the idea of comparing a numer- 



/ _q \ 1/4 / \ 1/2 ical solution with simple model nonlinear pen- 

Rmin = ( "j— 7 To~T ) 7^ • (5) 

\ 47r l9o + ^ 7r J/ V s '/ dulum, one may ask if it is possible to obtain 

From the expression for the effective poten- a better accord with the numerical experiments 
tial for i?o they obtained dependence of the using different ansatzes. We study this question 
monopole moment < r > and the breathing- in Section 3, and it seems that only the station- 
mode frequency on parameters gi, O. The ary Townes soliton can be fit accurately, but not 
frequency of small oscillations (breathing mode) the stabilized breather solutions, 
around the minimum is given by |4j 



8^ i n , „^2 tc\ 2. Modulated Townes soliton 



"br= vr(50 + 2vr) 2 . (6) 



Their numerical calculations were done for A method based on modulated Townes soli- 
go = —2tt. One can see that theoretical pre- ton used in Ref. (j, |j| should be mentioned. 
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The Townes soliton is a stationary solution to where n = 2, 3 is the dimension of the problem, 

the 2D NLS equation with constant nonlinearity In 2D, dr = 2irrdr, and in 3D dr = Amr 2 dr. 

g cr . In our notations \g cr \ « 1.8627T as 5.85. This For all t, we have I\ = 1. For the remaining 

solution is unstable: if \g\ is slightly increased I{ one can write down the dynamical equations 

or decreased, the solution will start to collapse of motion as 
or expand correspondingly. If the value of g is 

close to g cr , one may search for a solution of the . . n — 2 

h = h, h = h, h = g h + gh, (9) 

problem with fast oscillating g in the form of a n 

■ _ nvr2 n [°° SlVfdargV n _ x , 

modulated Townes soliton, as described in Refs. 5 ~~ g J Q g r q^. r ar{w) 

. A solution is sought in the form of m, , e , • , . . . 

° ine system of equations tor the momenta is not 

T , . r ,,,,1 _ r , ,•<? closed because of is, and one should make some 



r 2 ■ approximation in order to close it. In Ref. 

5 = a(t) + ^, a = a~\ (7) 

4a was assumed that 



it 



where Rt represents amplitude of the Townes I^r 2 

arg^ = ^-, (11) 

soliton. Then, starting from the approxima- 2 

tion ©, one can derive the evolution equation le - the phase factor is proportional to r 2 (so that 

for a(t) and so determine the dynamics of the a § ain [t is a PA "based method) and the coeffi- 

system. Note that the approach is also PA- cient of Proportionality is given by the ratio of 

based. It is inevitable if we are to use one- J 3 and I 2 . Then the system © poses dynamical 

parameter self-similar trial function in the form mva nants [8JJ. 

of ty(r,t)\ = Af(r/a,t). Qi = 2 (I 4 - g I 5 )I 2 - I/f, (12) 

Q 2 = 2/ 2 n/2 / 5 . (13) 



3. The method of moments 



With the help of these invariants, the system 
Another PA-based method we would like to becomes 



mention here is the method of moments 



■°- i 2 i ( ; 2 ^ 2 ^ +9 *V ( M) 



IT v / l j 1 y n /2 

introduces integral quantities I\, I2, h, ■■ as 2 \ 2 1% 



h = I \^\ 2 dr, h = / r 2 \^\ 2 dr, 



>3C 



Introducing X(t) = \J I2 it) one obtains 



\ dr dr 

1 f°° ( n n A The equation is analogous to that obtained 
h =5/ |VV| 2 + -<7WH 4 Ur, 

• ; ° ^ ' by other PA-based methods. One can investi- 

1 -1 r \^d 

5 4 jg ' ' r ' gate the obtained equation (fT3|) using various 



methods of nonlinear dynamics. The simplest of optical solitons) 14]. In Ref. j(| 



the solution 



Kapitza averaging method can be used again, is sought as an expansion in powers of I/O (in 
but of course it is better to use rigorous averag- our notation): 



1>{r,t) = A(r,T k )+n- 1 u 1 (A,C)+n- 2 u 2 (A,C)+...., 

(16) 

with < Uk >= 0, where < ... > stands for the 
average over the period of the rapid modulation, 
are the slow temporal variables (k = 
0, 1, 2, ...), while the fast time is £ = £lt. Then, 
for the first and second corrections the following 
formulas were obtained: 



ing technique since modern averaging methods 
are available 20] which have been extensively 
used already in plasma physics, hydrodynamics, 
classical mechanics . The authors of Ref. P] 
fulfilled rigorous analysis of model Eq. El us- 
ing results of Ref. It is important to have 
in mind that the relation between the exact dy- 
namics of the full system and that of the model 
1)15(1 of the method of moments remains unclear, 
therefore one cannot determine sufficient condi- u\ = —i\fMi— < >}\A\ 2 A, 

tions for stabilization, etc. In Ref. [8( it was no- ^ = [g(r)— < g\ >]dr 

Jo 

ticed that the correspondence between numerical r ir „., , l2 . . . 2 ,* . ..o ,m 

u 2 = [fi 2 - < ^2>}[2i\A\ 2 A t + iA 2 A* t + A(\A\ 2 A)] 

simulation of full 2D GP equation and dynam- i 

- \A\*A( -[(/.!-< Ml >) 2 -2M] + 
ics of the model system ()15f) is not good. As it z 

is seen from Fig. 3 of Ref. [g|, neither the fre- j? ^ 2 ^ 2 

quency of slow oscillations nor the position of the ^ 2 ~ Jo ^ 1 < Ml >)^ s i M— ^(</z x > < Hi 

minimum of the effective potential is predicted TT . , n r 

r Using these results, the following equation was 

correctly. Nevertheless, we found that in numer- , , . , r ,-, n , . c , , A , ^ \ j 

obtained tor the slowly varying held Ayr, 1q), de- 

ical stabilized solutions magnitudes of Qi and ■ , , , ? -2 

° ^ rived up to the order ot il . 

Q2 are often well-conserved, i.e. they oscillate ^ 2 

about some mean value (see Section 4). ~ = +1^1 ^ + ^ M (tt) 

-3|A| 4 AA + 2\A\ 2 &{\A\ 2 A) + A 2 &{\A\ 2 A*)]. 

(18) 



B. Direct averaging of the GP equation 



The above equation was represented in the quasi- 
Hamiltonian form 



1 + 6M ( ^ Vl 4 



n 



Ref. [a| also explores the Gaussian variational 
approximation. Beside that, a very promising 
method of directly averaging the GP equation 
was investigated. It is based on an analogous 
method used for the one-dimensional NLSE with 
periodically managed dispersion (in the context ~~ 2^ + ^Af— |V(|^4| ) ■ \^®) 



H q = I dV 



dA 5H q 
~dt ~ ~JA* : 

|V 4| 2 -2M (^] \A 
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However, some contribution was missed while 
deriving Eq. (|18j). Let us take into account the 
third correction u^(AX)'- 

4>(r,t) = A(r,T A )+0-V(AC)+tt~ 2 «2(AC) 

+ n- 3 u 3 (A,c) + .- (20) 



Then, up to terms of order Q~ 2 it changes noth- 
ing in r.h.s of Eq. (|TH|) (spatial part), but it adds 
to l.h.s. of Eq. (|18|) an undetermined term 
Q~ 2 dus/dC ■ This term has the same order Q~ 2 
as the terms from the second correction. So we 
do not get here a consistent equation for the slow 
field A because we do not have a closed set of 
equations for the second-order corrections (third 
order correction becomes second order correction 
after differentiating in time), and so the quasi- 
Hamiltonian (|19|) contains an undetermined er- 
ror of the second order in fi~ . The influence 
of the contribution is not very clear but require 
additional investigation. Nevertheless, formally 
the omitted terms have the same order as those 
responsible for the creation of the effective po- 
tential. Having in mind how many difficulties 
arise in averaging of systems of ordinary differen- 
tial equations |2Q], the rigorous direct averaging 
of the GP equation constitutes a very interesting 
and challenging open problem, since in principle 
it could reveal a true periodic solutions in such 
oscillating objects. 



III. VARIATIONAL APPROXIMATION 
WITH NON-GAUSSIAN ANSATZES 

Here we try to investigate the system more 
accurately using some non-Gaussian ansatzes 
and see if it is possible to get more accurate 
theoretical estimates. One may be interested in 
three dynamical quantities of the system: the 
value of critical nonlinearity g cr , slow frequency 
of breathing oscillations of the stabilized soli- 
ton LObr, and minimum of the effective potential 
Rmin about which the expectation value of the 
monopole moment < r > oscillates slowly. 

Table 1 summarizes results of variational pre- 
dictions for the critical nonlinearity g cr and fre- 
quency of small breathing oscillations using sev- 
eral different ansatzes. Note that the phase de- 
pendence of a one-parameter trial function is 
not important for calculating g cr . It is under- 
stood that if we choose a trial wavefunction 
with its amplitude in the form of \t/j(r,t)\ = 
Af[r/a(t)], then we need to use a phase factor 
with quadratic r— dependence in order for the 
ansatz to be self-consistent (i.e., the mass cur- 
rent generated by the changing parameter would 
be incorporated in the phase factor of an ansatz). 
On the other hand, since amplitude part of the 
trial function is just an approximation, one may 
try to use other forms of phase factor with the 
same functional form of the amplitude. 

When predicting the frequency of breathing 
oscillations from the corresponding effective po- 
tential, it is easy to obtain the result for small 
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amplitude linear breathing oscillations (given in g cr just because it is decaying too fast at large 

Table 1), but in actual stabilized solutions ampli- r. The supergaussian trial function provide a 

i 

tudes of breathing oscillations are not so small, better approximation, namely g cr = 7r2hi2ln2 

It is possible to take into account anhar- which corresponds to the supergaussian wave- 

monicity of breathing oscillations. As was men- function with ry = rj T = 21n2 < 2. Previously 

tioned earlier, all PA-based anzatzes produce the supergaussian ansatz was used to fit station- 

the nonlinear pendulum R + (a + bsmftt)/R 3 , ary solutions of some nonlinear problems includ- 

with a corresponding effective potential having mg NLS equation in the context of BECs |l9(]. 

Rmin = (-J^) 1 ^ 4 > = <J&Q\a/b\, where The superposition of two Gaussians in the form 

2 2 

u br is the frequency of the small amplitude A exp(-^)Cosh(7^2-) also enables one to ob- 
breathing oscillations (near the bottom of tain some improvement: g cr « 5.883. The Se- 
the effective potential). For larger breathing canth ansatz 

oscillations the (anharmonic) breathing fre- , A r-n/A n\ 2i 

^=^m exp[lSiR ' R)r] 



quency will be amplitude-dependent: cu%™ h 



works better, with only one parameter it over- 



j , D 2 _ D 2 comes the above-mentioned two-parameter trial 



functions. A very good approximation is pro- 
vided by the simplest ansatz among all consid- 
ered: 



with k = yjjfc^, where x x = R(, x 2 = R 2 

(Ri,R 2 being the turning points), x% is the 

third root of the equation h = ^ + .J^ 5 . The 

magnitudes of x±, x 2 , X3, h can be determined 

from numerically obtained breathing oscillations %h = — i — ( 1 + \ e xo I — - — h iS(R, R)r 2 \ 

V 3R^ V 2RJ y \ 2R ) 
(but results depend on the choice of a particular (21) 

anzatz). Even this improvement is not helpful, It fits the Townes soliton adequately both at the 

simply because the parabolic approximation is origin and asymptotically at infinite r ( a pre- 

not valid. exponential multiplier is not so important as the 

Finding g cr only might be considered as an exponential factor ). The pre-exponential factor 

approximation to the stationary Townes soli- is needed in order to fulfill the boundary condi- 

ton by a trial function so that the mass cur- tion in the origin lim r ^o < 00 • Note that 

rent term equals zero and that a phase factor in the supergaussian ansatz the former condition 

may be skipped from the calculations. It is is not fulfilled, otherwise (if one included it in 

known that the Townes soliton ifa = e lt RT(r,t) a similar way) the result would be better at the 

at large r has asymptotic behavior for its am- cost of more bulky calculations. The accuracy of 

plitude in the form Rt ~ eT r l\[r ■ So that the prediction implies that ansatz (|2*T|) provides a 

Gaussian ansatz is not very good for finding very good approximation to the Townes soliton 
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at fixed R, and could approximately represent sion for the phase factor. We also try the su- 
the modulated Townes soliton when R is time- pergaussian ansatz with fixed rj and with non- 
dependent and the phase factor with parabolic quadratic phase dependence (which is unfortu- 
r— dependence is used in accordance with the nately not self-consistent trial function) ip(r, t) = 
continuity condition. Aexp — \ a + ib ) rnT i where A, a, b, and 77 are all 
After obtaining estimates for g cr , one can use functions of time, parameter rj is fixed at the 
the above-mentioned ansatzes in order to find an value of its Townes soliton-like solution 77 = r] T = 
effective potential, its minimum and frequency 21n2. We find that such modification drasti- 
of the breathing oscillations of the monopole call y changes dynamical parameters of the sys- 
moment about this minimum in the same way tem - Still, the resulting model is the same clas- 
as it was done for the Gaussian ansatz. We sical nonlinear pendulum as in the Gaussian ap- 
checked the Sech ansatz and the supergaussian proximation, but with different parameters. The 
with quadratic phase dependence. In the su- rigorous way to employ the two-parameter su- 
pergaussian ansatz the parameter 77 was fixed pergaussian ansatz is to let rj be a dynamical 
at the value of its "Townes soliton-like" solu- variable and construct a phase factor fulfilling 
tion n = ifr = 21n2. In such a way the varia- continuity condition for the trial function. One 
tional approximation with supergaussian ansatz could then obtain the two-dimensional effective 
resembles method of modulated Townes soli- potential within the same Kapitza approach, 
ton. However, we find that such trial function As a useful test of applicability of the su- 
seriously underestimate minimum of the effec- pergaussian anzatz, we determine the critical 
tive potential (i.e. the mean value about which number of attractive BEC in the 3D parabolic 
the monopole moment oscillates). Nevertheless, trap studied in Ref. [l2|. Their numerical re- 
the result of the Gaussian ansatz is even worse suit was JV^. = 1258.5, while the gaussian ap- 
since for g = 2tt it gives the diverging expres- proximation yields N% = 1467.7. We found 
sion for R min and zero for frequency of slow the supergaussian prediction to be very accurate 
breathing oscillations uJbr, a s mentioned in Sec- N^ r G = 1236.1. 
tion 1 and J]. A natural idea for remedy is 

to use two-parameter trial functions to repro- jy NUMERICAL RESULTS 

duce the non-parabolic phase factor dependence 

on r. In the supergaussian ansatz it can be Numerical calculations reveal the fact that 
done by considering 77 as a dynamical (time- stabilized solutions do not have parabolic phase 
dependent) parameter. The problem is that it factors in contradiction to all the methods con- 
is difficult to obtain the self-consistent expres- sidered in Section 2 (except the method of direct 
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TABLE I: Variational predictions for the properties of stabilized solutions. 



Ansatz 


Amplitude part 
of the anzatz 


9cri 

analytical 
expression 


9cri 

approximate 
value 


(linear prediction, 

LObr = Kbr^/gi) 


Gaussian 
Supergaussian 

Secanth 
Exponential 


Aexp(-2^) 

ASech(^) 
A(l + &)exp(-&) 


2 TT 

1 

7r2hi2 l n 2 

27r ln22l n 2+i 
4ln2-i 

144^ 
— TT 


6.283 
5.919 
5.863 
5.875 


\/f(5o + 27r) 

/i"|21n2+l+ go (^y^T-)l 

\ / Q 4 In 2 - 1 



averaging) . The calculations were done using ex- 
plicit finite difference schemes. We use explicit 
finite differences of second and forth order for 
spatial derivatives and 4-th order Runge-Kutta 
method for time propagation. We use meshes 
varying from 2000 to 10000 points, timesteps 
At = 0.0001 ~ 0.0004, and spatial steps Ar = 
0.02 ~ 0.04. In addition, we found that it is very 
important to use absorbing (imaginary) poten- 
tial at the edge of the mesh, in accordance with 
the conclusions of Ref. Without such an 

adsorbing potential, a wave reflected from the 
edge sometimes destroys the otherwise stable so- 
lution. 

Following P|, initially we start with a Gaus- 
sian wavepacket in a parabolic trap. Then the 
trap was slowly turned off while the oscillating 
nonlinearity was slowly turned on in a way sim- 
ilar to Ref. J]. In Figure ^one can see indeed 
the creation of a stabilized soliton. In Figure 
\& and \& oscillations of amplitude of the wave- 
function at the origin are shown. It decays very 
slowly. In fact, this is in accord with the calcu- 
lations of Ref. 41: after a careful examination 



of the corresponding figures in that paper one 
notices the same behavior. Monopole moment 
grows very slowly (Figs. ^,,c). We checked that 
in the case when the trap is not turned off com- 
pletely, the norm is conserved during the same 
long time with a high accuracy (of order 10 -8 ), 
so decay is certainly not due to numerical errors. 

In Fig. [21 configuration of the quasi-stabilized 
wavefunction is shown. One can see the smooth 
core pulse profile, tiny oscillations in the tail, 
and an outgoing cylindrical wave leaking from 
the core pulse. In Fig. |2^ the behavior of the 
phase factor is shown. It is seen to differ from 
parabolic with r considerably. 

Figs. |2J,g shows the slow decay of the norm 
of the solution due to the flux of atoms from 
the core to infinity. We made a series of nu- 
merical experiments with different parameters. 
We found that the behavior of the matter-wave 
pulse is often unpredictable. When the Gaussian 
approximation predicts stabilization, in the cor- 
responding numerical solution it does not nec- 
essarily occur. Neither can the method of mo- 
ments give reliable predictions for the stabiliza- 
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tion. We checked the latter method carefully. 
As it was mentioned already in Sections 1 and 2, 
the method relies on the crucial approximation 
of Eq. (fTTl ) . It is due to this approximation one 
obtains the existence of dynamical invariants Q\ 
and Q2 (see Eq. H5)l . As a result, dynamics is de- 
termined by Eq. (|15[). Returning back to Figure 
12 we see a snapshot of the phase factor, arg tp, 
of a stabilized solution. It clearly demonstrates 
that none of the PA-based methods reproduce 
the dynamics of the system adequately. Only 
at small r the parabolic law is fulfilled, while 
the deviation from this quadratic dependence is 
very strong even at r < 1, where the amplitude 
of the solution is not small at all (and is sufficient 
to drastically influence the dynamics of the sys- 
tem). Snapshots at other moments produce sim- 
ilar results: the phase of the solution is chang- 
ing with time but remains very far from being 
parabolic in r. It is easy to check that dynami- 
cal properties of the system within a variational 
approximation are very sensitive to r— depen- 
dence in the phase factor of a trial function. To 
check the dynamics further, we calculated time 
evolution of the "invariants" Q\ and Q2 in the 
stabilized solution. They are constants in the 
model but not in the exact numerical solution. 
We found that in the numerical quasi-stabilized 
solution these magnitudes oscillate around some 
mean value. Actually, it was already found in 
Ref. j3] that the method of moments does not 
work for Gaussian initial data, still it is inter- 
esting to trace dynamics of relevant magnitudes. 



The time evolution of Q\ and Q2, and other mag- 
nitudes related to the method of moments are 
shown in Figures and El It is seen that the 
magnitudes of Q\ and Q2 related to a stabilized 
soliton undergo slow oscillations. 

When calculating values of Q\ and Q2, and 
other properties of the quasi-stabilized solution 
it is necessary to stop integration at some reason- 
able value of r = r max (we take r max = 20 where 
the amplitude of the wavefunction becomes very 
small (of order 10~ 4 in our case). In that way 
we separate the properties of the quasi-stabilized 
soliton from that of the tail which, although has 
very small amplitude, can carry large moments 
I2 , ^3 and would give large contribution to Q\ 
and Q2 (so that in the corresponding figures 
we presented these quantities for the core soli- 
ton and the whole solution (including tail) sep- 
arately) . 

Similar features can be seen in Fig. where 
calculations with g$ = —7.0 are presented. Sev- 
eral snapshots of the phase factor at different 
moments are presented in order to demonstrate 
that non-quadratic behavior of the phase factor 
is typical. Time evolution on very long time is 
traced. We find that sometimes magnitudes of 
Qi and Q2 of stabilized solutions are almost con- 
served (undergoing small oscillations about its 
mean value) despite the strongly non-quadratic 
behavior of the phase factor. It suggests that the 
method of moments developed in 8] might pro- 
vide useful perspective for studying the problem 
and it would be fruitful to extend it taking into 
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account non-parabolicity of the phase factor. part in nonlinear optics. As was studied in Ref. 

3], in the periodically alternating Kerr media 



V. CONCLUDING REMARKS 



the stabilization of beams is possible. Mathe- 
matically, one deals with a similar NLSE. In- 
Despite there are many publications dedi- stead of the time-dependence of the scattering 
cated to the stabilization of a trapless BEC length of BEC one has dependence of the media 
by the rapidly oscillating scattering length, it nonlinearity coefficient on the coordinate z along 
seems that the strong non-parabolic behavior which a beam propagates: 
of the phase of the stabilized wavefunction has 

not been brought to attention yet. It should be i 

iu z + -V 2 r u + ~/{z)\u\ 2 u = 0, (22) 
noted that the role of deviation of the phase pro- ^ 

file of NLSE solutions from the parabolic shape where the diffraction operator V 2 r acts on 
was addressed previously in thecontexts of soli- the transverse coordinate x and y. Nonlinear- 
tons in optical fibers in Refs. [17j, [l^] ■ ity coefficient 7(2) jumps between constant val- 



Despite that several independent methods ues j± of opposite signs inside the layers of 

were used previously, we have seen that three widths L±. The analysis of this problem was 

of the four theoretical methods used rely on the done using variational approximation based on 

unwarranted parabolic approximation, while the a natural Sech ansatz U = A(z) exp[ib(z)r 2 + 

fourth method (direct averaging of GP equation) i(f)(z)]Sech[r/w(z)]. However, behavior of the 

is, strictly speaking, incorrect, despite its inspir- phase factor was not checked aposteriori . We 

ing motivation (in the sense that the omitted see that it would be useful to investigate the 

terms has the same order as those responsible for problem of (2+l)-dimensional solitons in a lay- 

the creation of the effective potential ). Besides, ered medium with sign-alternating Kerr nonlin- 

we find that there is no evidence presently for earity in a greater detail because behavior of the 

stabilization in a strict sense. It seems that the phase factor of the numerical solution has not 

numerical examples presented so far deal with been reported yet. Interplay between the phe- 

quasi-stable solutions which slowly decays due to nomenon of stabilization in Kerr media and BEC 

the leaking of atoms from the core pulse as an was addressed also in Ref. 22| in the context of 

outgoing cylindrical wave. It means that even stabilization of (3+1)- dimensional optical soli- 

from a numerical point of view there are no evi- tons and BEC in periodic optical-lattice poten- 

dence for true 2D solitons (breathers) yet. tial (without addressing the issue of validity of 

It should be mentioned also that the phe- the parabolic approximation), 

nomenon of BEC stabilization has its counter- Returning back to the BEC stabilization, we 
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note that the two main difficulties should be re- 
solved in the future: the non-trivial behavior of 
the argument of the stabilized wavefunction, and 
the possibility to stop the leak of atoms from the 
tail of the solution. 

Using several non-Gaussian variational func- 
tions, we were able to determine accurately one 
of the magnitudes characterizing the stabiliza- 
tion phenomena: critical nonlinearity g cr , but 
not other dynamical properties such as the fre- 
quency of slow oscillations. 



VI. ACKNOWLEDGEMENTS 

Alexander Itin (A.I.) acknowledges support 
by the JSPS fellowship P04315. A.I. thanks 
Professor Masahito Ueda and Professor S.V. 
Dmitriev for helpful discussions. 

This work was supported in part by Grants- 
in-Aid for Scientific Research No. 15540381 and 
16-04315 from the Ministry of Education, Cul- 
ture, Sports, Science and Technology, Japan. 



[1] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. 

Stringari, Rev. Mod. Phys. 71, 463 (1999). 
[2] I. Towers, B. A. Malomed, J. Opt. Soc. Am. B 

19, 537 (2002). 
[3] L.V. Hau, M. M. Burns, J. A. Golovchenko 

Phys. Rev. A 45, 64686478 (1992). 
[4] H. Saito and M. Ueda, Phys. Rev. Lett. 90, 

040403 (2003). 
[5] F. Kh. Abdullaev, J.C. Bronski, and R. M. 

Galimzyanov, cond-mat/0205464; Physica D 

184,319 (2003). 
[6] F. Kh. Abdullaev, J.G. Caputo, R. A. Kraenkel, 

B.A. Malomed, Phys.Rev. A 67, 013605 (2003). 
[7] F. Kh. Abdullaev, A.Gammal, L.Tomio, T. 

Frederico, Phys. Rev. A 63, 043604. 
[8] G. D. Montesinos, V.M. Perez-Garcia, P.J. Tor- 
res, Physica D 191, 193210 (2004). 
[9] Gaspar D. Montesinos et. al, Chaos 15, 033501 

(2005). 

[10] S.K. Adhikari, Phys. Rev. A 69, 063613 (2004). 
[11] M. Matuszewski et. al, Phys. Rev. Lett. 95, 

050403 (2005). 
[12] C. Huepe, S. Metens, G.Dcwcl, P. Borckmans, 



M.E. Brachet, Phys. Rev. Lett. 82, 1616. 
[13] Y. Kagan, EL. Surkov and G.V. Shlyapnikov, 

Phys. Rev. A 54, R1753 (1996) 
[14] T.S. Yang and W. L. Kath, Optica Letters 22, 

985(1997). 

[15] J. Lei, M. Zhang, Lett. Math. Phys. 60, 9 
(2002). 

[16] V.M. Perez-Garcia, H. Michinel, J. I. Cirac, M. 

Lewenstein, and P. Zoller, Phys. Rev. Lett. 77, 

5320 (1996); Phys. Rev. A 56, 1424 (1997). 
[17] Opt.Commun. 147, 317 (1998). 
[18] Progr. Optics 43, 71 (2002). 
[19] D. Anderson, M. Lisak, A. Bertson, Pramana 

J. Phys. 57, 917 (2001). 
[20] V.I.Arnold, V.V.Kozlov, and A.I.Neishtadt, 

Mathematical aspects of classical and celestial 

mechanics (second edition, Encyclopaedia 

of mathematical sciences 3, Springer- Ver lag, 

Berlin, 1993). 
[21] Sec, for example: A.P.Itin, A.I.Neishtadt, 

A.A.Vasiliev, Phys. Lett. A. 291, 133 (2001); 

A. P. Itin, R. de la Llave, A. I. Neishtadt, A. 

A. Vasiliev, Chaos 12, 1043 (2002); AP. Itin, 



15 



A. A. Vasiliev, A.I. Neishtadt, Physica D 141, [22] S. K. Adhikari, Phys. Rev. E 71, 016611 (2005). 
281 (2000). 



16 



<r> 



<r> 




0.90 J , , , , , , . ,— 

100 104 108 112 116 



time 

(e) 

FIG. 1: (a) Oscillations of the monopole moment after turning off the trap. Parameters are g = —2n, 
gi = 8n, ft = 30. The trap was turned off completely at T Q ff — 30. (b) Time evolution of amplitude 
of wavefunction at the origin, (c) Oscillations of the monopole moment on longer time scale, (d) Time 
evolution of amplitude of wavefunction at the origin on longer times, (e) The oscillations of the monopole 
moment from previous figure on finer scale. Tiny high frequency oscillations are seen. 
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FIG. 2: Configuration of the quasi-stabilized wavefunction. Parameters are the same as in previous figures, 
(a) A snapshot of an amplitude profile, (b) Tiny oscillations in the tail of the quasistabilized solution (c) 
Amplitude of the wavefunction far from the origin (the tail plus outgoing cylindrical wave), (d) Real part 
of the wavefunction far from the origin multiplied by -Jr. (e) Snap-shot of the phase factor of the quasi- 
stabilized solution. It can be seen that it is parabolic only at very small r. The curve has an inflection point 
at r < 1. f),g) The slowly decaying norm of the solution. Although the trap was turned off at t — T a ff = 30, 
the norm remains almost constant until the flux of atoms leaking from the core soliton reach the edge of the 
mesh and begin to disappear. After that it decreases slowly. 
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FIG. 3: Oscillations of the monopole moment, (a) go — —6.5, ft = 35, g\ = 10n . Initial frequency of the 
parabolic trap is u(0) = 0.8. (b) g = —6.5, = 30, g\ = 14.5 . Initial frequency of the parabolic trap is 
w(0) = 1. Quasistabilizcd solution is destroyed after several oscillations. 
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(c) (d) 

FIG. 4: Time evolution of the integral quantities Qi,Q2- (a) Oscillations of "shortened" Q 2 (designated 
as 92)- We integrate expressions entering Eq. from r = to r « 8 so that it characterizes the core 

part of the solution (quasi-stabilized soliton) without the oscillating tail, (b) Time evolution of full Q2- 
The expressions (|13|l were integrated from r = to r w 120 so that it includes large contribution from the 
oscillating tail, (c) Time evolution of "shortened" Q x (designated as q-J. We integrate expressions entering 
Eq. from r = to r w 8 so that it characterizes the core part of the solution. Dynamics of the core 
soliton for quite a long time is almost independent of the behavior of the tail which after reaching the edge 
of the grid begin to disappear, (d) Time evolution of the full Q\ (including large contribution from the 
oscillating tail which depends on location of the absorbing potential and the mesh size ). 
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FIG. 5: Time evolution of the moments (l2,h,Ii)- (a) Oscillations of the second moment < r 2 > of the core 
soliton (designated as 12). The boundary of the core of the quasi-stabilized soliton was taken to be r w 8. 
(b) Time evolution of the second moment I 2 =< r 2 > of the whole solution including tail (this magnitude 
depends on mesh size, here v max w 120) (c) Time evolution of I4. (d) Time evolution of I3. 
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FIG. 6: Oscillations of the monopole moment of a quasistabilized solution. Parameters are go = —7.0, 
il = 40, #1 = 8tt . Initial frequency of the parabolic trap was chosen to be w(0) = 4.0. (a) Time evolution of 
the monopole moment on very long time, (b) Detailed picture of the time evolution of the monopole moment 
of a quasistabilized solution about t 600. (c) Detailed picture of the time evolution at t = 1800 <~ 2000. (d) 
Decaying norm of the solution, (e) Time evolution of the integral quantity Qi (calculated for the core part of 
the wavefunction) (f) Time evolution of Qi. (g)Several snap-shots of the phase factor of the quasi-stabilized 
solution (made at different moments). Note that typical behavior of the phase factor is not quadratic with 
r at all. 



